PFU counting
Installation
Install packages if not installed yet
# Install package. Enter "a" in console when prompted
if (!require("BiocManager", character.only = TRUE)) install.packages("BiocManager")
if (!require("EBImage", character.only = TRUE)) BiocManager::install("EBImage")
if (!require("tidyverse", character.only = TRUE)) install.packages("tidyverse")
if (!require("broom", character.only = TRUE)) install.packages("broom")Load packages
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✓ ggplot2 3.3.5 ✓ purrr 0.3.4
## ✓ tibble 3.1.6 ✓ dplyr 1.0.7
## ✓ tidyr 1.1.4 ✓ stringr 1.4.0
## ✓ readr 2.1.1 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(broom)
library(EBImage)##
## Attaching package: 'EBImage'## The following object is masked from 'package:purrr':
##
## transposeNote
Be careful that this Rmd and the accompanied Rscript only work for the cropped image (see below) that only show ONE community at ONE dilution factor.
This code only works for counting the total number of colony. I am not able to distinguish the different colony morphology yet. But I think for those with >100 colony this should help speed up a bit.
Image processing
Read a test image. I have cropped the data such that
# Read an image
img <- readImage(here::here("data/raw/test2.jpeg"))
display(img, method = "raster", all = T)
Reduce the dimension of the image from 3 to 1. Using the most contrast dimension
# Grey scale
colorMode(img) = Grayscale
display(img, method = "raster", all = T)
img <- img[,,3]
display(img, method = "raster")
Blur the image as the pre-treatment for gating out the background noise
# Blur
img_bl <- gblur(img, sigma = 1)
display(img_bl, method = "raster")
Set the global threshold
# Global threshold
threshold <- otsu(img_bl)
img_th <- combine(mapply(function(frame, th) frame > th, getFrames(img), threshold, SIMPLIFY=FALSE))
display(img_th, method = "raster")
Label the objects. Each object is a group of adjacent pixels.
img_label = bwlabel(img_th)Gating out the background. Be cautious that here I am setting a manual threshold (e.g., objects number of pixel <= 10 are considered noise).
# Remove background
temp <- img_label %>% table()
index_rm <- names(temp)[which(temp <= 10)] # Pixel size for gating out background
img_clean <- rmObjects(img_label, index_rm)
display(img_clean, method = "raster")
Watershed separates the adjacent colonies
# Watershed
img_ws <- watershed(distmap(img_clean), 1)
display(colorLabels(img_ws), method = "raster")
You can get shape data for each colony: area, perimeter, radius, etc. The number of objects are the number of colonies.
colony_shape <- computeFeatures.shape(img_ws) %>% as_tibble()Number of total colonies
nrow(colony_shape)## [1] 34colony_shape## # A tibble: 34 × 6
## s.area s.perimeter s.radius.mean s.radius.sd s.radius.min s.radius.max
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 533 106 12.8 2.00 8.84 17.3
## 2 420 103 11.3 1.74 8.13 15.5
## 3 290 74 9.15 1.20 6.30 11.8
## 4 272 90 8.86 1.83 4.54 12.2
## 5 258 76 8.77 1.30 6.29 11.4
## 6 248 66 8.39 1.04 6.41 10.8
## 7 243 68 8.40 1.06 6.09 10.6
## 8 270 81 8.86 1.39 5.88 11.9
## 9 332 94 10.2 2.60 4.60 14.9
## 10 228 72 8.09 1.25 4.44 10.5
## # … with 24 more rowsTest on another image
# Read an image
img <- readImage(here::here("data/raw/test3.jpeg"))
display(img, method = "raster", all = T)
# Grey scale
colorMode(img) = Grayscale
display(img, method = "raster", all = T)
img <- img[,,3]
display(img, method = "raster")
# Blur
img_bl <- gblur(img, sigma = 1)
display(img_bl, method = "raster")
# Global threshold
threshold <- otsu(img_bl)
img_th <- combine(mapply(function(frame, th) frame > th, getFrames(img), threshold, SIMPLIFY=FALSE))
display(img_th, method = "raster")
img_label = bwlabel(img_th)
display(img_label, method = "raster")
# Remove background
temp <- img_label %>% table()
index_rm <- names(temp)[which(temp <= 10)] # Pixel size for gating out background
img_clean <- rmObjects(img_label, index_rm)
display(img_clean, method = "raster")
# Watershed
img_ws <- watershed(distmap(img_clean), 1)
display(colorLabels(img_ws), method = "raster")
In this case, the difference in colony size is very large. It makes sense to tell them apart by the size
colony_shape <- computeFeatures.shape(img_ws) %>% as_tibble()
colony_shape## # A tibble: 50 × 6
## s.area s.perimeter s.radius.mean s.radius.sd s.radius.min s.radius.max
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 6069 368 44.1 3.73 37.4 53.5
## 2 5467 317 41.4 2.79 35.4 49.4
## 3 3517 242 33.6 4.76 24.8 41.9
## 4 2616 218 28.8 4.07 18.5 38.7
## 5 2664 254 28.9 1.31 26.0 32.5
## 6 948 190 17.1 5.63 3.47 27.3
## 7 261 66 8.70 0.961 7.10 11.4
## 8 237 59 8.27 0.909 6.26 10.8
## 9 236 60 8.27 0.813 6.19 10.2
## 10 241 58 8.39 0.746 6.73 9.73
## # … with 40 more rowsRelative count by clustering
colony_shape %>%
as.matrix() %>%
kmeans(2) %>%
broom::tidy() %>%
pull(size)## [1] 45 5